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Abstract 

In this paper, we propose an edge detection technique based on some local smoothing 
of the image followed by a statistical hypothesis testing on the gradient. An edge 
point being defined as a zero-crossing of the Laplacian, it is said to be a significant 
edge point if the gradient at this point is larger than a threshold s(e) defined by: 
if the image / is pure noise, then P(||V/|| > s(e) | AI = 0) < e. In other words, a 
significant edge is an edge which has a very low probability to be there because of 
noise. We will show that the threshold s(e) can be explicitly computed in the case of 
a stationary Gaussian noise. In images we are interested in, which are obtained by 
tomographic reconstruction from a radiograph, this method fails since the Gaussian 
noise is not stationary anymore. But in this case again, we will be able to give 
the law of the gradient conditionally on the zero-crossing of the Laplacian, and 
thus compute the threshold s(e). We will end this paper with some experiments 
and compare the results with the ones obtained with some other methods of edge 
detection. 
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1 Introduction 



This work is part of some specific physical experiments which consist in study- 
ing radially symmetric objects [3]. These objects are composed of several ma- 
terials and one of the interesting features is the location of the frontier between 
the different materials. 



To describe such an object, it is enough to give the densities of the materials 
on a slice of the object that contains the symmetry axis. An example of studied 
object is given on Figure 1. 



Fig. 1. Slice of a studied object. 



To look at the interior of this object, a radiography is performed (see Figure 
2(a)), then a tomographic reconstruction is computed (Figure 2(b)) and finally 
an edge detection is made (Figure 2(c)). 
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(a) (b) (c) 

Fig. 2. (a) Radiograph of the object of Figure 1. (b) Tomographic reconstruction, 
(c) Edge detection on the tomography. 



At this point, let us mention that the tomographic reconstruction is not an 
usual one. Indeed, the usual inverse Radon transform (and the usual recon- 
struction algorithms such as filtered back-projection) operates on a slice of the 
object that is orthogonal to the rotation axis. Here, because of the radial sym- 
metry assumption, the reconstruction can be global [5]. This reconstruction 
will be detailed in Section 4.1. 

As we can see on Figure 2(c), many detected edges do not correspond to 
real features. This is due to the high level of noise. For the time being, the 
selection of the edges is manually executed. The goal of this work is to perform 
this selection automatically. For that purpose, the edge detector will not be 
changed but we will compute also other significant features that will allow us 
to select the "true" edges. 

The ideas used here mainly come from previous work of Desolneux, Moisan and 
Morel [4]. Informally speaking, they define the notion of significant edges by 
computing the probability that some edge-related events appear in an image 
of pure noise. When this probability is small enough, the edge is probably a 
feature of the image and not due to the noise. Unfortunately, their method 
assumes that the noise is stationary which, as easily seen on Image 2(b), is 
not the case in our study because of the tomographic inversion (see Section 
4.5.1 for some examples of results obtained with their method). Moreover, 
their study is quite general and, apart from the stationarity, no assumption is 
made on the noise. 

In our case, as we deal with specific images, the noise is well-known and some 
statistical models can be used. Indeed, we may suppose that the noise on the 
radiograph (2(a)) is a Gaussian white noise with mean zero and with a variance 
that can easily be estimated. Then, a tomographic inversion is performed. As 
this operation is linear, we still obtain a Gaussian noise but it is now correlated 
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and non-stationary. 



The edge detector will not be modified here. It consists in estimating the 
Laplacian at each point, and edge points are then defined as the zero-crossings 
of the Laplacian. As we already said, we only add some features that will 
discriminate the significant edges. The easiest feature to compute is a kind of 
contrast measurement C based on a gradient estimate. To be more precise, we 
consider an image I of pure noise (that is a realization of our model of noise 
after tomographic reconstruction), estimate the gradient and the Laplacian of 
/ at a point (u, v) (with an abuse of notation, we will denote by V/(ii, v) and 
AI{u,v) these estimates and by C(u,v) the contrast value) and we compute, 
for a fixed e > 0, the smallest value s(e) for which 

P(C(u, v) > s(e) | AI(u, v) = 0)< e. (1) 

Then, we perform an edge detection on the studied image / (where we also 
estimate V/ and Af by the same method) and we keep the points (u, v) of 
the studied image / that satisfy 

• Af(u,v) = (an edge is present at point (u,v)). 

• C(u,v) > s(e) (this edge is significant). 

/From a statistical point of view, this consists in performing an hypothesis 
test. We consider a point (u,v) where an edge takes place (A(u,v) = 0) and 
we test the null hypothesis "the zero-crossing of the Laplacian is due to the 
noise" . The level of the test e is arbitrarily chosen and related to the number 
of false detections allowed. It will be set to e — 1CT 5 hereafter. Let us mention 
that the threshold value s(e) varies slowly with respect to e. For instance, in 
the case of a white noise (see Section 3), the threshold value can be computed 
explicitly and is proportional to \J — lne. When the null hypothesis is rejected, 
the edge is retained as it comes from a "true" feature of the image, whereas 
when the null hypothesis is accepted, the zero-crossing of the Laplacian may 
come from the noise and the edge is not meaningful. 

Let us mention at this point that such statistical approaches have already been 
used for edge detection in [11], [10] or [8]. They usually use estimates of the 
gradient based on finite differences which fail in our case. Moreover, the noise 
is in most cases stationary. Let us also cite [2] where the authors have modified 
the method of [4] to take into account the non-stationarity of some images, 
by a local noise estimate. Their work is still general and does not make any 
assumption on the noise structure. As we deal with specific experiments, the 
noise is always the same and well-known and we can take proper advantage 
of this knowledge. 

The paper is organized as follows: in Section 2, we present the edge detector 
based on the estimate of the gradient and the Laplacian. Then, in Section 
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3, our method is presented in the case of a Gaussian white noise. Of course, 
this does not correspond to our case but the computations are easier and 
show the performance of this method. In Section 4, we will first describe the 
tomographic inversion and the operators involved, and then describe the noise 
model we have to deal with. We will then apply the significant edges detection 
method in the framework of this non-stationary noise. We will end the section 
with some experiments and comparisons with other methods. 



2 Estimating the Gradient and the Laplacian 

In this section, we introduce a method for edge detection. We consider that 
the image is a real- valued function (u,v) h- > f(u,v) of two continuous real 
parameters u and v. Then, we say that there exists an edge at point (u, v) if the 
Laplacian of / is zero at this point. Moreover, the computation of the contrast 
function C will be based on the gradient of / (see the end of this section for the 
choice of this function). As the images are very noisy, these derivatives cannot 
be estimated by usual finite differences. The method used here, sometimes 
known as Savitsky-Golay smoothing, consists in locally approximating the 
image by a polynomial. The derivatives of the polynomial are then identified 
with those of the image. 



2.1 An optimization problem 

Let (u, v) denote the point where we want to compute the first and second 
order derivatives of the image /. We choose 2 parameters : d which is the 
maximum degree of the approximating polynomial and r which is the radius 
of the ball on which we perform the approximation. We denote by B r (u,v) 
the ball of radius r centered at point (u,v). We will simply write B r when 
the center of the ball is the origin (0, 0) of 1R 2 . We are then looking for a 
polynomial P of degree less that d such that 

E(P) = j B (f(u + x,v + y)-P(x,y)) 2 dxdy (2) 

is minimal among all polynomials of degree less than d. In other words, we are 
looking for the best approximation of / by a polynomial of degree less than d 
on the ball B r (u, v) in the sense of the L 2 -norm. 

This is an optimization problem where the unknowns are the coefficients of 
the polynomial. As the problem is convex, there is a unique solution (given 
by the orthogonal projection of / on the space of polynomials of degree less 
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than d) which is easily computed by solving the equations 




where the Oj's denote the coefficients of the polynomial. 

Role of the ball radius. Two parameters are arbitrary chosen in this method. 
The first one is the ball radius r. The larger r is, the more effective the smooth- 
ing is. The influence of the noise is therefore attenuated with a large r but 
the location of the edge is then less precise. We must consequently make a 
balance between noise smoothing and edge detection accuracy. For instance, 
if we have a small level of noise or if the edges are very complicated (with high 
curvature), we must choose a small value for r. 

Role of the polynomial degree. The second parameter is the polynomial degree. 
Here again a large value of d gives a better approximation but does not smooth 
the noise enough. In fact, as we are, in a first step, interested in the points 
where the Laplacian is zero, it appears that a second-order polynomial is 
enough. Of course, the estimate of the first order derivatives with a polynomial 
of degree 2 is not very good and highly depends on the size of the window B r . 
But we will see that this drawback can be useful for the choice of a contrast 
function. 

In what follows, the approximation is made with a polynomial of degree d = 2, 
and the first and second order derivatives of the image are identified with those 
of the approximating polynomial. 

2.2 Computations with a second order polynomial 

Let us first introduce some notations. In the following, we will set 



As the ball B r is symmetric, we have that by(r) = as soon as i or j is odd 
and that bij(r) = bji(r) for all In order to have simple expressions, we will 
also set: 



Lemma 1 The gradient and the Laplacian of the polynomial of degree 2 which 
is the best approximation of f on the ball B r (u,v) for the L 2 -norm, being 
respectively denoted by V r f(u,v) = (^-(u,v),^-(u,v)) and A r f(u,v), are 






and 
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given by: 



drf ( v 
9 r f , \ 

A r f(u,v) 



77~t I x f(u + x,v + y)dxdy 

b(r) J B r 

b(rj Ib y f( u + x > v + y) dxd y 



P(r) Jb.. 



f(u — x, v — y) (a(r) + x 2 + y 2 ^j dx dy. 



Proof : 

We consider a polynomial of degree 2 which we write 

P(x, y) = a xx x 2 + a yy y 2 + a xy xy + a x x + a y y + a . 

The equations obtained by writing VE(P) = 0, where E{P) is given by Equa- 
tion (2), are: 



ho(r)a xx + b 2 2(r)a yy + b 20 (r)a = x 2 f(u + x, v + y) dxdy 

J B r 

b 22 (r)a xx + b 40 (r)a yy + b 20 (r)a = y 2 f(u + x,v + y) dxdy 

J By 

b 22 (r)a xy = / xy f(u + x,y + v) dxdy 

J B r 

b 20 (r)a x = / xf(u + x,y + v)dxdy 

b 20 (r)a xx + b 20 (r)a yy + b 00 (r)a = / f(u + x,v + y)dxdy 

J B r 



We then obtain the following estimates for the derivatives: 



(0, 0) = a x = , 1 / x f(u + x,v + y)dxdy 
b 20 {r) Jb t 

(0,0) = a y = —— y f(u + x,v + y)dxdy 
bon(r) JB r 



dP 
dx 
dP 
dy 

AP(0,0) = 2(a : 



J2Q 
.).)■ T" Q>yy) 



ho + b 22 - 2 -^ 



[ f(u + x,v + y) (— + x 2 + y 2 ) dxdy 

JB r V b 00 I 



□ 
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2. 3 Choice of the contrast function 



We would like to use a contrast function based on the estimates of the first 
and second derivatives of the image / obtained in the previous section. 

The simplest contrast function we can choose is the norm of the gradient: 

C 1 (u,v) = \\V r f(u,v)\\. 

Indeed, the value of this norm tells how sharp the edge is. This contrast 
function is efficient and will be used when the images we deal with are piecewise 
constant. 

However, in many cases, the objects we handle are not homogeneous and their 
images contain some slopes (see Figure 3). In this case, the gradient norm is 
not a good contrast function. Indeed, let us consider an image with a constant 
slope with some noise (see Figure 4). We would like to say that no edge is 
significant in that case. However, the value of the gradient norm (which will 
be close to the value of the slope) will always be greater that the threshold 
value s when the noise level is small. 




Fig. 3. Object with an inhomogeneous material 

In the latter case, we take advantage of the dependence of the first order 
derivatives estimates with respect to the ball radius. Indeed, the estimates of 
the gradient in the case of the constant slope in Figure 4 will not depend on 
the size of the window (see Figure 4) whereas, when an edge (a discontinuity) 
occurs, the estimates do depend on that radius (see Figure 5). So, we can use 
as a contrast function the function 

C 2 {u,v) = \\V r J{u,v) -V r J{u,v)\\ 

where r\ < r 2 and V r / denotes the value of the gradient estimate with a ball 
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of radius r. 




Fig. 4. A noisy constant slope: the gradient of the approximating polynomial does 
not depend on the radius r. 




Fig. 5. An edge on a noisy slope: the approximating polynomials with two different 
values of the radius r. 



3 Significant edges in the case of a Gaussian white noise 

3.1 White noise and Wiener integral 

We recall here the definition and the main properties of a white noise in a 
continuous setting and of the Wiener integral. We refer to [12], [6] or [7] for 
more on white noise and the Wiener integral. 

Definition 1 A Gaussian white noise on M 2 of variance a 2 is a random func- 
tion W defined on the Borel sets A ofM, 2 of finite Lebesgue measure (denoted 
by \A\) such that 

• is a Gaussian random variable (r.v.) with mean and variance o~ 2 \A\ 

• If Ax fl A 2 — 0, the r.v. W(A 1 ) and W(A 2 ) are independent and 

W(A X U A 2 ) = W{A X ) + W{A 2 ). 
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Such a function W exists but is not a true measure since the two-parameters 
process 

B(s,t) := W((0,s] x (0,tf 



(usually called the Brownian sheet) is of unbounded total variation. 

Nevertheless we can define the so-called Wiener integral / fdW for every func- 
tion / in L 2 (R^_). We can also define the derivatives of the Brownian sheet in 
the sense of Schwartz distributions (although the Brownian sheet is nowhere 
differentiable). Thus, we define 

and we have 

/ fdW = / f(u,v)B(u,v)dudv a.s. 
for every function / in the Schwartz space. 

With a slight abuse of notations, we call B a Gaussian white noise and we 
always denote by J R 2 f (u, v) B(u, v)du dv the Wiener integral with respect to 

this white noise, for every function / e L 2 . The main properties of this integral 
are 

• For every /, the r.v. / f(u, v)B(u, v)du dv is a Gaussian r.v. with mean 

• For every /, g, the random vector 



f(u,v)B(u,v)dudv, / g(u,v)B(u,v)dudv 

is Gaussian with covariance 

a 2 J ^ f(u, v)g(u, v)du dv. 

We will use these properties to compute the laws of V/ and AI 



3.2 Laws of the gradient and of the Laplacian 



We suppose here that our noise is a Gaussian white noise, of variance a 2 . As 
we have already said, this case is not the one we are interested in and our 
method is probably over-performed by other standard methods in that case. 
The goal of this section is to present our method in a simple case where the 
computations are easy to do and can be carried out in a continuous setting. 
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We will only focus here on the case of piecewise constant objects and therefore 
we will use the contrast function C\. 



Lemma 2 // the image I is a Gaussian white noise of variance a 2 , then 



dx ' dy ' r , 



is a Gaussian vector with mean zero and covariance matrix 



b(r) 

b(r) 








^ V(r,a) ) 



0~ 2 2 

where V(r, a) = / (a(r) + x 2 + y 2 ) dxdy. 

p (r) J B r v 



Proof : We compute the laws of the approximate derivatives of / when I = B. 
We recall that these derivatives are given by 

d r I . . 1 r ■ . 

= ITT / xB(u + x,v + y)dxdy 
ox b(r) J B r 



I »■<'.) - 1J^\J B yB(u + x,v + y)dxdy 



d r I 

dy ~' " y b(r) Jb t 
1 r 

A r I(u, v) = / .B(-u + x, w + y) ( a(r) + :r 2 + y 2 ) cte dy. 



(3(r) Jb 

Because of the stationarity of B, they have the same law as 
d r I 1 



(0, 0) = — — x B(x, y) dx dy 



dx b(r) Jb t 

8 r I 1 



a 1 I f 

7T~( ' ) = ITT / yB{x,y)dxdy 
dy b(r) JB r 

A r /(0, 0) = J-f B(x, y) (a(r) + x 2 + y 2 ) dx dy. 

p[r) JB r v ' 



As we deal with Wiener integrals, we deduce that the vector 



d A Q A A T 

s dx dy j 



is a Gaussian vector with mean zero. 
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To compute its covariance matrix, let us recall that, if X and Y are random 
variables defined by 



then we have 



X= h 1 (x,y)B(x,y)dxdy 

J B r 

Y= h 2 (x,y)B(x,y) dxdy 

J B r 



Cov(X, Y) = a 2 / hi(x,y)h 2 (x,y) dx dy. 

J Br 



Consequently, we have for instance: 



Covf^,M 



a 



dx ' dy J b 2 (r 



■ xydxdy = 0. 

J B r 



By some analogous calculations, we finally get the following covariance matrix 
for our Gaussian vector: 

'sty N 

b(r) 



T7~T 

b(r) 

[ V(r,a) ) 



where 



2 

V(r, a) = -^^y jT (a(r) + x 2 + y 2 ) rfxrfy. 



Thanks to this lemma, we immediately have the following properties: 



□ 



The random variable ||V r /|| is the sum of two squared independent Gaus- 
sian random variables which have the same variance. It is therefore dis- 
tributed as a x 2 Taw. More precisely, its law is 



b(r) 



X 2 (2) 



where x 2 (2) denotes a x 2 -law with two degrees of freedom. 

The random variable A r J is a Gaussian random variable with mean zero 

and variance V(r, a). 

The random variables ||V r /|| 2 and A r I are independent. 
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3. 3 Computation of the threshold 



Proposition 1 Let I be a Gaussian white noise and let s(e) be the threshold 
value such that 

p(||V r /|| > s{e) | A r I = 0) < e. 

Then 

2a 2 



Ins. 



Proof : To begin with, as the random variables || V r /|| 2 and A r I are indepen- 
dent, we can forget the conditioning and only compute 

F( llVrJll > s(e)) = Pf II Vrlf > s(s) 2 



As a consequence of Lemma 2, we have that the law of ||V r /|| 2 is ^yx 2 (2)- 
Now, since the density of a x 2 (2) law is the one of a T (|, l) law, we have that 
the law of ||V r /|| 2 is given by 



F(l|V,J|| 2 > S 2 )=/^ 2 ie-^t = exp (- 



b(r)s 2 



2a 2 J ' 

This finally leads to the announced threshold value s(e). □ 
3.4 Experiments 



We consider the piecewise constant object of Figure 1 with some additive 
Gaussian white noise. The densities of the different materials of this object 
are: 

• 1 for the outer material, 

• 0.8 for the material inside the circle, 

• 0.3 for the other inner material. 

The standard deviation of the Gaussian noise is a — 0.2 in the experiments 
of Figure 6 and is a = 0.4 in the experiments of Figure 7. Both images have 
the same size 512 x 512 pixels. The experiments have been carried out with a 
ball of radius r = 12 pixels. 

The different images of Figures 6 and 7 are respectively: 



(a) The noisy image. 
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• (b) The zero-crossings of the Laplacian with the contrast function G\ visual- 
ized in grey-level (the white color corresponds to high value for the contrast 
function C\). 

• (c) The extracted significant edges (e = 10~ 5 ). 




(a) (b) (c) 

Fig. 6. (a) The noisy image {a = 0.2). (b) The zero-crossings of the Laplacian with 
the contrast function C\ visualized in grey-level, (c) The extracted significant edges 
(e = 10" 5 ). 




(a) (b) (c) 

Fig. 7. (a) The noisy image (a = 0.4). (b) The zero-crossings of the Laplacian with 
the contrast function C\ visualized in grey-level, (c) The extracted significant edges 
(e = 10" 5 ). 



In the case of a signal-to-noise ratio large enough (Figure 6), all the edges are 
well detected and the "false" edges are removed. Let us nevertheless mention 
that, with our method, the edges which have a high curvature are smoothed. 
This drawback is even more important when the ball radius r is large (the 
influence of the value of this radius will be studied in the experiments of the 
next section). 
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When the noise level is rather large (Figure 7), some edges of the image cannot 
be extracted from the noise (it happens when the contrast associated to this 
edge is close to the noise level). 



4 Significant edges in the case of a Gaussian white noise on the 
radiograph 

4-1 Tomography 

Let us turn now to the more realistic case we are interested in. As we mentioned 
it in the introduction, we first make a radiography of an object. Tomography is 
the inverse problem associated with reconstructing the initial object from its 
radiograph. This is now a well-known problem as it is the key tool in medical 
scanner imagery (or other medical imaging systems). 

To begin with, let us describe what a radiography is from a mathematical point 
of view. The studied object is exposed to X-rays that go through it. Some of 
the X-photons are absorbed. As an output, we observe the quantity of X- 
photons that have not been absorbed by the material, and we thus measure 
in some sense the "mass" of material the ray went through. More precisely, 
if the object is described by its density fi (which is a function of the space 
coordinates), what can be measured at some point of the receptor is 



where "ray" means the straight line that goes from the source to the studied 
point of the receptor (we suppose that the X-rays source is just a point, which 
implies that the previous line is unique). 

We also assume that the X-rays source is far away from the object so that 
the rays are assumed to be parallel. Then, to reconstruct any object from its 
radiographs, we must turn around the object and make a radiography for every 
angle 9 G [0, 7r). This leads to the so-called Radon transform of the object, 
which is known to be invertible. This is the principle of the medical scanner. 

In our case, as the object is radially symmetric, if we turn around the object 
with for rotation axis the symmetry axis of the object, all the radiographs 
are exactly the same. Consequently, a single radiograph of such an object is 
enough to perform the tomographic reconstruction. Indeed, if f(x, y) denotes 
the density along a slice that contains the symmetry axis (see Figures 1 and 




ray 
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8), then a radiograph of this object is given by 



g(u, v) = 2 f(x, v) dx. 

J\u\ \/X z — U 



This is a linear transform and we will denote it hereafter by 

9 = Hf. 

As we already said, this linear operator H is invertible and we in fact know 
explicitly its inverse on the space of continuously differentiable functions g: 

1 f+°° 1 dg 

f(x,y) = {H g)( x ,y) = / - = — (u,y)du. 

vr Jx v« z — x z ou 




Fig. 8. Radiography of a radially symmetric object. 



Our assumption on the noise is that it is an additive Gaussian white noise on 
the radiograph (i.e. on g). But what we want is to study the object given by 
/. So we must transform the white noise by the operator if -1 . Unfortunately, 
because of the singularity of the integral at i = 0, we cannot apply the 
operator H^ 1 to a white noise B, even in a L 2 -sense. Therefore, we will work 
in a discrete framework: the images / and g are naturally discretized (as they 
are numerical images). This leads to a discretization of the operator H, which 
we will still denote by H and which now may be viewed as a matrix. The 
discretization is made in such a way that the symmetry axis (x = 0) is settled 
between two pixels so that the previous singularity does not appear. This 
matrix is then invertible and we denote by H^ 1 its inverse which we can make 
now operate on a discrete Gaussian white noise. 
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4-2 Law of the noise on the tomographic reconstruction 



Let us consequently consider a field r] = (Vi,j)i<i<p,i<j<n of i.i.d. random Gaus- 
sian variables with mean and variance a 2 . Let us define I = (Iij) the random 
field obtained after tomographic reconstruction i.e. after making H~ x operate 
on rj = (rjij). In fact, as the X-rays are supposed to be parallel, the recon- 
struction can be made line by line independently and therefore, if we consider 
the row vectors 

ffi = (?7i,i,...,?7i,n) and U = . . . , I i>n ) 

then, there exists an invertible matrix M (independent of i, and of size n x n) 
such that 

h = ffiM. 

Consequently, the law of / is characterized by the following properties: 

• I = (Iij) is a Gaussian random field. 

• For i ^ k, Ii and 1^ are independent. 

• For each i, the vector Jj is a Gaussian vector of mean and covariance 
matrix 

r = ct-m 1 m. 

where M* denotes the transpose of M. 



4-3 Laws of the gradient and of the Laplacian 



The expressions obtained in Section 2 for the gradient and for the Laplacian of 
an image in a continuous setting are easily translated in the discrete framework 
we now deal with. Indeed, we have 

d r I 1 

v ' (%,j)eB r 
d r I ( , 1 -r 

A r I(u, v) = ( a ( r ) + {2 + i 2 ) tu+i^+j 

P\ r > (i,j)eB r 

where B r now denotes the discrete ball of radius r i.e. 

B r = {(i,j), * 2 +f<r 2 } 

and where the constants a(r), f3(r), b(r), . . . are the discrete analogous of the 
constants of Section 2. 
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With these estimates, the contrast functions C\ and C 2 are easily comptuted. 
They are both of the form 

C(u,v) = y/C*(u,v) + C*(u,v) 

with 

C x (u,v) = jcjjl u +i,v+i and C y (u,v) = ^ iCjjI u+ i tV+ j, 

i,j i,3 

where the coefficients are given by: 

(1) In the case of the contrast function Ci, 

Cij = b(r) 1(i ' j)eBr ' 

(2) In the case of the contrast function C 2 with two balls of radius r 1 < r 2 , 

° ij = b(^j 1{i ' j)eBri ~ b(^j 1{l ' j)eBr2 ' 

Therefore, the computations of the laws will be similar and they will be treated 
simultanously using the coefficients c^-. 

When the contrast function C2 is used with two radii r\ < r 2 , we then compute 
the Laplacian A r I with the larger ball radius, that is with r = r 2 . 

Lemma 3 For both contrast functions C\ and C 2 , the vector 

(C x (u, v),C y (u, v),A r I(u, v)) 

is a Gaussian vector with mean and covariance matrix of the form: 

o-l 

In particular, we have that C y is independent of (C x , A r I) . 

Proof : The lemma is a consequence of the two following remarks. The first 
one is that, in both cases for the contrast function, the coefficients are 
symmetric: = c_jj and = Cj ; _j. Thus they satisfy that whenever k or I 
is odd then 

E <Vcy = (3) 

(i,j)€B r 
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The second remark is that the vectors Ii and Ik are independent if i 7^ k. And 
we thus have 

'0 ifi^k, 



E [Ii,jh,i] - 



r(j,i) xi = k. 



We can now compute the covariance matrix. For instance, let us start with: 
E [C X Cy] = E jkcijCkM[I u +i,v+jIu+k,v+i} 

(i,j,k,l) 

= E 3iCijCuT{v + j,v + I) 

= E J r ( v + i. t, + 0E ic iO C H = °- 

(3,1) i 



Similar computations give E [C y A r I] = and 



a 2 :=E 



a 2 A :=E 



Cl 



(Kiy 



E jlc ij c i iT(v + j,v + I); 
E i 2 CijC U r(v + j,v + I); 

1 



<r xA :=E[W] = 



E (a(r)+i 2 +j 2 )(a(r) + < 2 + Z 2 )r(t; + j J T; + 0; 
P \ r ) (i,j,i)en r 

E JC^a(r) + 2 2 + Z 2 )r> + j,t; + Z), 



^2) (ij.ijen^ 



where we have set Q r = I) such that G £> r and (i, /) G S r }. 



□ 



4-4 Computation of the threshold 



Now, as we have no more independence between the first and the second 
order derivatives we must compute the conditional law of the contrast function 
knowing that A r I = 0. 

Proposition 2 Let C be one of the two contrast functions. Then, the random 
variable \\C\\ 2 is distributed, conditionally on {A r I = 0}, as the sum of the 
square of two independent Gaussian random variables, with mean zero and 
respective variance 



a 2 and a 2 , A=0 



2 : 
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that is a Gamma law with parameters \ and \{o~ 2 + °"x|A=o)- 

The threshold value s(e) defined by 

P(||C|| > s(e) | A r / = 0) <e 

can no longer be computed explicitly but a numerical approximation is easy 
to get as the Gamma density is well-known. 

Proof : C y is independent of the pair (C x , A r I). Thus, conditionally on {A r I = 
0}, the random variables C y and C x are still independent and the conditional 
law of C y is the Gaussian distribution with mean and variance a 2 . 

Now, if D 2 := o 2 x a\ — a 2 A 7^ 0, then the law of the pair (C x ,A r I) has a 
density which is given by 

/x,A(*i,*2) = 2^e-3^'^^ t 
where A is the inverse of the covariance matrix, i.e. 

1 



A =D 2 



( a \ 



\ -<T X ,A °l ) 



Let us recall that, if /a denotes the Gaussian density of A r I, then the law of 
C x conditionally on A r I = has a density given by 

/s,A(*l,0) 

/a(0) 

and so is Gaussian with mean zero and variance 

D 2 

a x\A=0 — 9 ' 
a A 

This result is still valid when D = since it implies that C x and A r I are 
proportional and thus the law of C x conditionally on A r I = is Gaussian 
with mean and variance (it is not random anymore). □ 



4-5 Experiments 



4-5.1 Case of a piecewise constant object 

To begin with, we still study the piecewise constant object of Figure 1 de- 
scribed in Section 3.4. Let us recall that this image represents a slice of the 
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object that contains the symmetry axis. The 3-dimensional object is obtained 
by rotation around the vertical axis that goes through the middle of the image. 

In that case, we will use the contrast function C±, which is simply the norm of 
the gradient. The experiments of Figure 9 correspond to a ball radius r = 12 
pixels. 

We start with the image of the radiograph obtained after the application of 
matrix H to our initial image. Then a Gaussian white noise is added to this 
radiograph. Then tomographic inversion (application of the matrix H^ 1 ) is 
performed. This gives the image of Figure 9(a). As we already mentioned it, 
the noise is not stationary, it is now correlated and its variance depends on 
the distance from the symmetry axis. For instance, if the standard deviation 
of the Gaussian white noise on the radiograph is a = 4, the variance of the 
noise on the tomography is about 2a 2 = 32 near the axis, 0.02<r 2 = 0.32 at a 
distance of 65 pixels from the axis and 8.1CTV 2 = 0.128 at the edge of the 
image located on the right at 200 pixels from the axis. 




(a) (b) (c) 

Fig. 9. (a) Reconstructed object from a single noisy radiograph, (b) Contrast value 
at the zero-crossings of the Laplacian for the contrast function C\, (c) Significant 
edges. 

We notice that the edges are not significant near the symmetry axis; the noise 
is too important here in order to extract the true edges from the noise. Let us 
add that the smaller the difference of the densities of the material is, the larger 
the region where the edges are not significant around the axis is. Even when 
the edges are significant, the noise and the method used to detect them can 
lead to noisy edges. Moreover, some details are lost because of the smoothing 
due to the size of the ball. 

Let us compare the results obtained with different ball radii (see Figure 10). 
When the ball radius is small, the edges are more accurate but some are not 
significant: the smoothing of the noise is not enough to get rid of it. On the 
contrary, when the radius is large, most of the edges are detected but small 
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details are lost because of this smoothing. 




r = 6 r = 12 r = 20 

Fig. 10. Significant edges obtained with different ball radius: from left to right: 
r = 6, r = 12 and r = 20 

Since the edges separate two materials, one included in another, they must be 
closed curves. Usually, an operator has to close them manually. Our method 
gives open edges. It does not mean that there is no edge between the materials: 
it simply means that the noise level is too high to give an accurate position 
of the edge. Therefore, we can then close the curves manually, or by usual 
curve completion methods, but this will not tell which closure is better (i.e. 
the closest to the real shape). 

Comparison with other methods. We will give here the results obtain with two 
other methods which have both the advantage of directly providing closed 
curves. 

• The first method is the one introduced in [4] . One keeps only the meaningful 
level lines of the image, which are defined by: the minimum of the norm of the 
gradient along the level line is larger than a threshold T(e). This threshold 
is computed from the gradient histogram of the image. The meaning of this 
definition is that such curves have a probability less than e to appear in a 
pure noise image (with same gradient histogram as the original image). The 
results obtained with this method are shown on Figure 11. On the first row: we 
smooth the image of Figure 9(a) by convolution with a Gaussian kernel with 
respective standard deviation 2 and 4 pixels. And then, on the second row, 
we have the respective obtained meaningul level lines. This experiment clearly 
shows that, since the noise model is not adapted to the image (in particular, the 
non-stationarity is not taken into account), many false contours are detected. 

• The second method is the famous Mumford-Shah segmentation for piecewise 
constant images [9]. Given an observed image g defined on a domain D, one 
looks for the piecewise constant approximation g of go that minimizes the 
functional 

Ex(g)= J D \g - g \ 2 + XLength(K(g)), 
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where Length \K(g) j is the one-dimensional measure of the discontinuity set of 
g (which is a set of curves denoted by K(g)) and A is a parameter which weights 
the second term of the functional. The results obtained with this method are 
shown on Figure 12 for three different values of A. The main drawbacks of 
this method are: (a) there is no verification that the obtained contours are not 
due to the noise; (b) the parameter A has to be fixed, and the results are very 
dependent on its value. 





Fig. 11. First row: the image of Figure 9(a) is smoothed by convolution with a 
Gaussian kernel with respective standard deviation 2 and 4 pixels. Second row: the 
meaningul level lines of each image. 
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Fig. 12. Results obtained with the Mumford-Shah segmentation for piecewise con- 
stant images, for three different values of A. From left to right, the number of regions 
in the segmented image is respectively 3, 6 and 7. 



4-5.2 Case of an inhomogeneous material 



Let us turn now to a more realistic case: the materials are not homogeneous 
and consequently the object is no more piecewise constant (see Figure 13). 
As already said, the use of the contrast function C\ fails in that case. This is 
illustrated by Figure 14. In this image, one can notice that there are many false 
detections especially in the parts of the image where it is not constant. Figure 
15 gives the significant edges obtained with the contrast function C2 with two 
ball radii T\ — 6 and r2 = 12. With this contrast function, we eventually get 
only the "true" edges. 




(a) (b) (c) 

Fig. 13. (a) Inhomogeneous object, (b) Its noisy radiograph, (c) Tomographic re- 
construction 
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Fig. 14. Significant edges with the C\ contrast function: there are many false detec- 
tions. 




Fig. 15. Significant edges with the C2 contrast function: only the "true" edges are 
obtained. 
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